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A  numerical  method  was  developed  for  calculating  the  Rayleigh  Surface  Wave  (RSW) 
velocity  in  arbitrarily  oriented  single  crystals  in  360  degrees  of  propagation.  This  method 
relies  on  the  results  from  modern  analysis  of  RSW  behavior  with  the  Stroh  formalism 
to  restrict  the  domain  in  which  to  search  for  velocities  by  first  calculating  the  limiting 
velocity.  This  extension  of  existing  numerical  methods  also  leads  to  a  natural  way  of 
determining  both  the  existence  of  the  RSW  as  well  as  the  possibility  of  encountering  a 
pseudo-surface  wave.  Furthermore,  the  algorithm  is  applied  to  the  calculation  of  elastic 
properties  from  measurement  of  the  surface  wave  velocity  in  multiple  different  directions 
on  a  single  crystal  sample.  The  algorithm  was  tested  with  crystal  symmetries  and  single 
crystal  elastic  moduli  from  literature.  It  was  found  to  be  very  robust  and  efficient  in 
calculating  RSW  velocity  curves  in  all  cases. 
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1.  Introduction 

Surface  acoustic  waves  (SAW)  such  as  Rayleigh  surface  waves  (RSW)  are  important  in  a  number  of  industries  and  appli¬ 
cations  such  as  nondestructive  evaluation,  materials  characterization  and  acoustic  microscopy.  These  waves  propagate  with 
distinct  particle  motions  due  to  the  mixing  of  both  shear  and  compressional  modes,  and  thus  also  have  a  different  velocity. 
Furthermore,  in  an  anisotropic  material,  the  wave  velocity  and  particle  displacements  change  as  a  function  of  propagation 
direction.  This  is  often  represented  with  a  slowness  curve  which  gives  the  reciprocal  of  velocity  as  a  function  of  propagation 
direction. 

One  application  area  for  SAW’s  is  in  material  property  characterization.  It  has  been  shown  in  the  past  that  measurement 
of  RSW  velocity  on  a  single  crystal  In  multiple  directions  of  propagation  can  help  calculate  the  single  crystal  elastic  constants 
of  a  material  [1],  More  recently,  modern  experimental  techniques  have  been  used  to  measure  the  velocity  of  the  surface 
waves  much  quicker  than  previously  possible  by  relying  on  laser  excitation  and  detection  of  elastic  strains.  In  [2]  the 
RSW  velocity  measurements  were  used  to  characterize  the  orientation  map  of  a  sample  by  calculating  the  velocity  given 
known  orientations  and  elastic  constants  and  then  inverting  the  response  to  get  to  orientation.  In  [3],  the  response  over 
the  entire  surface  of  the  sample  was  used  to  characterize  the  elastic  properties  of  a  material  by  generating  a  distribution 
of  velocities  in  the  many  different  orientations  present  in  a  polycrystalline  aggregate.  In  all  of  these  studies,  a  key  step  is 
accurately  calculating  the  RSW  velocity  given  known  orientation,  propagation  direction,  and  elastic  constants.  These  works 
all  employed  the  methods  of  Farnell  [4]  for  numerical  computation  of  the  RSW  velocity. 
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While  numerical  methods  have  been  used  extensively  in  the  past  for  calculation  of  the  RSW  velocity,  there  are  a  number 
of  authors  who  have  derived  various  equations  for  analytically  calculating  RSW  velocities  [5,6],  These  derivations  often 
rely  on  the  Stroh  formalism  for  elastodynamic  equations,  although  not  always  [7],  This  formalism  reorganizes  the  stiffness 
matrix  in  an  elastic  wave  problem  and,  after  application  of  boundary  conditions,  results  in  an  eigenvalue  problem  that 
combines  the  governing  differential  equation  with  the  boundary  condition  to  determine  the  wave  propagation  constant, 
ultimately  leading  to  a  polynomial  equation  whose  roots  are  directly  related  to  the  wave  velocity.  With  an  application 
of  rotation  matrices  to  the  elastic  stiffness  matrices,  this  problem  could,  in  theory,  be  solved  at  any  arbitrary  orientation 
and  propagation  direction.  However,  the  eigenvalue  problem  that  results  from  the  decomposition  of  the  stiffness  matrix  is 
not  analytically  tractable  in  an  arbitrarily  rotated  crystal  also  having  arbitrary  crystal  symmetry.  Even  so,  authors  have  used 
other  methods  and  linear  algebra  techniques  to  argue  past  the  eigenvalue  problem  and  arrive  at  a  new  determinant  equation 
for  calculating  the  secular  equation  for  RSW  velocity  explicitly  [8],  This  was  done  for  an  arbitrary  elasticity  tensor  in  [9], 
in  which  a  determinant  equation  that  was  a  function  of  only  the  elastic  constants  and  wave  velocity  was  used  to  derive 
various  analytical  secular  equations  that  had  previously  been  derived.  Though  this  approach  gives  a  very  general  method 
of  calculating  the  coefficients  of  the  secular  equation  explicitly  as  a  function  of  elastic  constants  and  rotation  angles,  the 
generalized  version  of  the  equations  are  very  difficult  to  derive,  even  in  symbolic  math  languages  such  as  Mathematica.  As  a 
result,  the  equations  are  either  solved  for  specific  orientations,  propagation  directions,  and  crystal  symmetries  or  numerical 
techniques  are  used  to  solve  the  equations  [10], 

While  these  analytical  computations  are  not  tractable  for  a  generalized  elasticity  matrix,  there  are  aspects  of  the  deriva¬ 
tion  that  can  make  the  numerical  solutions  more  stable  and  quicker.  In  this  paper,  a  numerical  technique  for  solving  for 
the  RSW  velocity  and  pseudo-surface  wave  (PSW)  velocity,  as  well  as  determining  the  existence  of  the  RSW  in  an  arbi¬ 
trary  direction  on  an  arbitrarily  oriented  crystal  is  presented.  The  method  is  based  on  the  typical  iterative  search  procedure 
described  in  [4].  In  that  approach,  the  velocity  space  is  searched  over  with  non-gradient  based  optimization  methods.  In 
a  situation  in  which  the  solver  has  found  real  eigenvalues  (leading  to  non-decaying  waves  into  the  solid  half-space,  or 
non-surface  wave  behavior),  the  real  eigenvalue  that  gives  the  lowest  value  of  the  boundary  condition  determinant  is  used. 
This  leads  to  many  different  minima  that  become  difficult  to  sort  through  in  arbitrary  cases.  It  also  leads  to  discontinuous 
gradients  of  the  boundary  condition  determinant  with  respect  to  velocity,  preventing  the  use  of  gradient-based  optimization 
routines.  The  typical  approach  to  solving  this  problem  is  to  perform  the  inverse  many  times  and  determine  which  minima 
corresponds  to  the  RSW  velocity.  This  results  in  a  non-trivial  solution  time,  so  much  so  that  the  forward  model  is  not 
used  to  directly  compute  the  elastic  constants,  and  a  database  approach  is  typically  employed  [2,11],  An  alternate  method 
for  calculating  the  velocities  is  to  perform  simulations  of  the  entire  wave  field  using  numerical  simulations,  e.g.  [12-15], 
These  methods  give  a  much  more  in  depth  view  of  the  physics  of  the  waves  interacting  with  a  complex  anisotropic  material 
and  can  be  used  to  model  the  incident  wave  field,  but  they  require  significant  computational  resources  relative  to  simply 
calculating  the  velocity  using  FarnelTs  method. 

In  the  new  numerical  method  described  in  the  following  sections,  the  limiting  velocity  is  first  found  using  an  iterative 
search  procedure.  This  provides  a  solid  lower  and  upper  limit  to  the  search,  which  allows  the  use  of  quickly  converging 
search  methods  to  find  the  RSW  velocity.  Furthermore,  by  comparing  the  RSW  velocity  with  the  limit  velocity,  conditions 
under  which  pseudo-surface  waves  exist  and  dominate  the  response  of  the  material  can  be  determined.  The  method  was 
tested  on  many  different  crystal  symmetries  and  directions  of  propagation,  and  the  advantages  and  disadvantages  are  dis¬ 
cussed.  An  alternative  method  of  finding  the  RSW  velocity  is  then  discussed  that  is  based  on  finding  the  surface  impedance 
matrix.  This  method  is  slower  but  is  much  more  robust  in  actually  converging  to  the  RSW  velocity  when  local  minima  are 
obscuring  the  answer  in  the  boundary  condition  determinant. 

2.  Elastic  wave  propagation  at  the  boundary 


For  the  sake  of  consistent  notation,  the  fundamentals  of  wave  propagation  in  elastic  media  will  be  given  here.  The  goal 
of  the  analysis  is  to  determine  the  surface  wave  velocity  in  an  elastic  media  given  a  certain  propagation  direction  along  a 
surface  cut  along  an  arbitrary  crystallographic  plane.  Thus,  a  form  of  a  solution  that  would  indicate  surface  wave  propagation 
at  an  unknown  velocity  is  substituted  in  to  the  governing  equations  of  elasticity  and  traction-free  boundary  conditions  are 
then  applied,  and  the  velocity  that  satisfies  the  boundary  conditions  is  considered  a  true  surface  wave  velocity. 

The  governing  equations  of  motion  in  a  linear,  anisotropic,  homogeneous  medium  with  no  force  terms  are  stated  using 
summation  notation  as: 


3t2  “  9xjX/ 


In  this  equation,  p  is  the  mass  density,  Uj  is  the  displacement  in  the  ith  direction,  and  Cijid  is  the  elastic  stiffness  tensor. 
The  solution  to  this  equation  that  would  coincide  with  a  wave  propagating  on  the  surface  of  the  sample  can  be  written  as: 


_  gg-ik(m-x-l-pn-x-vt)  ^2) 

Here,  a  is  polarization  vector  that  describes  the  particle  displacements  of  the  wave,  i  =  m  is  a  unit  vector  in  the 

direction  of  propagation  of  the  wave,  n  is  an  outwardly  pointed  unit  vector  normal  to  the  surface  of  the  sample  on  which 
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the  wave  propagates  and  is  orthogonal  to  m,  p  is  a  complex  constant  to  be  determined,  v  is  the  velocity  of  the  wave,  and 
k  is  the  wave  number.  Note  that  any  point  inside  the  half-space  results  in  n  •  x  being  negative.  This  equation  is  simply 
a  statement  that  there  is  a  wave  with  velocity  v  traveling  in  the  m  direction  on  the  surface  of  a  half-space  that  has  the 
outward  normal  vector  n.  Since  m  is  specified  as  a  real  unit  vector,  it  travels  infinitely  in  the  m  direction  with  no  decay. 
However,  the  addition  of  complex  p  in  front  of  the  n  ■  x  term  indicates  there  is  both  an  oscillating  term  as  well  as  a  decay 
term  (assuming  lm(p}  <  0)  in  this  direction,  which  is  crucial  for  a  surface  wave  to  exist.  As  stated  earlier,  the  surface  wave 
cannot  propagate  into  the  sample,  thus  It  must  decay  with  respect  to  the  distance  away  from  the  surface.  Furthermore,  the 
thickness  of  this  layer,  and  thus  the  value  of  the  imaginary  component  of  p,  is  dependent  on  the  wave  velocity.  However, 
the  wave  velocity  is  unknown  at  this  point,  which  means  p  and  v  must  be  determined  simultaneously. 

One  way  to  solve  this  equation  would  be  to  plug  (2)  Into  (1).  Using  the  procedure  shown  in  [5],  we  can  define  the  3x3 
matrices: 


Q={Q/,k} 


3 

i,J=l 


R={Ri,k} 

'r={Ti,k] 


3 

^  Cijkimjnt 

!,i=i 

3 

Cywnjn, 


(3) 


In  these  expressions,  Si  k  is  a  Kronecker  delta  and  m,  and  n,  is  are  the  ith  component  of  the  m  and  n  vectors,  respectively. 
These  can  be  used  to  simplify  the  resulting  expression  from  plugging  (2)  into  (1}  down  to: 


[q-Fp(r-FR^) -Fp2T]a  =  0  (4) 

where  a  is  the  vector  of  constants.  This  is  the  quadratic  eigenvalue  problem  that  results  from  satisfying  the  equations  of 
elasticity  for  a  plane  wave,  with  p  being  the  eigenvalue  and  a  the  eigenvector.  At  this  point,  if  the  wave  velocity  was  known, 
it  could  be  used  to  solve  this  eigenvalue  problem  numerically  to  calculate  the  a,  and  p  and  thus  know  the  wave  profile 
entirely. 

The  Rayleigh  wave  velocity  is  still  unknown  at  this  point.  A  further  step  of  satisfying  the  boundary  conditions  must  be 
taken  to  calculate  it.  The  traction  at  the  surface  n  •  x  =  0  is  expressed  as: 


t  =  -ikbe-‘‘'<"’-*-''f> 
where  the  vector  b  is  given  by 


(5) 


b  =  (r^  -F  pt)  a 


(6) 


Equations  (4)  and  (6)  can  be  combined  to  give  the  new  6x6  eigenvalue  problem  known  as  the  Stroh  eigenvalue  problem: 


where 


(7) 


'^=  _(2+RT  Ir’^  _rT  1 

This  is  a  linear  eigenvalue  problem  for  calculating  the  eigenvector  a  and  the  eigenvalue  p.  However,  there  are  6  eigenvalues 
and  6  eigenvectors  that  satisfy  this  equation.  In  the  range  of  velocities  between  0  and  the  limiting  velocity  (discussed  in 
Section  3.1),  the  eigenvalues  are  all  complex  and  come  In  conjugate  pairs.  If  the  3  eigenvalues  and  eigenvectors  satisfying 
lm{p}  <  0  are  labeled  with  the  index  a  =  1, 2,  3,  the  traction  at  the  surface  (5)  can  then  be  written  as  a  linear  combination 
of  the  three  eigenvalue/eigenvector  pair  solutions: 


3 

t  =  _ike-ik(moc-vt)  ^  0 

a=l 

where  Ca  are  arbitrary  complex  constants  that  are  not  all  equal  to  zero.  This  can  be  rewritten  as: 


[bi 


^2 


b3] 


Cl 

C2 

C3 


=  Bc  =  0 


(9) 
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This  implies  that  the  boundary  condition  at  the  free  surface  is  satisfied  when  the  matrix  of  eigenvectors,  B  is  singular, 
ultimately  leading  to  the  expression: 


det{B}  =  0 


(10) 


The  velocity  that  satisfies  this  condition  satisfies  both  the  governing  equations  as  well  as  the  boundary  conditions,  and 
thus  constitutes  a  viable  surface  wave  velocity.  In  a  numerical  scheme  the  velocity  is  iterated  over  until  the  condition 
(10)  is  satisfied.  Algorithms  of  this  type  have  been  suggested  since  as  early  as  [16],  though  even  earlier  references  allude 
to  computational  frameworks  to  calculate  the  Rayleigh  wave  velocity  numerically.  However,  with  the  surge  of  research  in 
applying  the  mathematical  rigor  of  the  Stroh  formalism  to  the  calculation  of  Rayleigh  wave  velocities  in  the  2000’s,  new 
results  are  available  that  make  the  numerical  algorithms  much  more  stable  and  quick,  which  is  absolutely  critical  for  using 
them  in  the  context  of  materials  characterization. 

Another  important  concept  in  the  analysis  of  surface  waves  is  the  surface  impedance  matrix,  defined  as: 


(11) 


Z=-iBA 


where  A  is  the  matrix  formed  by  combining  the  displacement  vectors,  a  into  a  single  matrix.  The  inverse  is  guaranteed 
to  exist  in  the  event  that  (4)  has  three  distinct  eigenvalue  solutions  and  complex  conjugate  solutions.  In  the  event  that 
this  equation  results  in  repeated  eigenvalues,  generalized  eigenvectors  can  be  computed,  and  these  vectors  are  used  in 
forming  the  matrix  A.  In  either  case,  the  eigenvectors  are  independent,  and  thus  the  matrix  inverse  exists.  The  surface 
impedance  matrix  will  be  used  when  direct  calculation  of  RSW  velocities  with  (10)  is  not  feasible.  In  the  next  sections, 
a  numerical  algorithm  is  detailed  which  is  very  stable  for  calculation  of  the  entire  curve  for  Rayleigh  waves  on  the  surface 
of  an  arbitrarily  oriented  single  crystal. 

3.  Numerical  scheme 

Rayleigh  surface  wave  velocities  can  be  calculated  by  finding  the  zeros  of  equation  (10)  above.  The  algorithm  starts  at  a 
certain  velocity,  assembles  the  sub-matrices  (3),  calculates  the  eigenvectors  and  eigenvalues  from  (7),  selects  the  eigenvec¬ 
tors  that  predict  a  decaying  solution  (i.e.  negative  imaginary  components  in  the  eigenvalues)  and  assembles  tbe  matrix  B 
and  calculate  its  determinant.  The  velocity  is  then  iteratively  changed  until  the  zero  is  found. 

3.1.  Limiting  velocity  calculation 

There  are  some  difficulties  that  arise  when  trying  to  take  advantage  of  the  modern  computation  routines  for  1-D  root 
finding.  For  instance,  there  is  only  a  finite  region  where  the  eigenvalues  of  N  are  complex  and  come  in  conjugate  pairs.  Even 
building  the  B  matrix  outside  of  this  range  is  not  straight  forward  for  a  computational  routine.  The  typical  solution  to  this 
problem  is  that  found  in  Farnell  [4],  where  the  real  eigenvalue  that  gives  the  lowest  absolute  value  of  the  determinant  of  B 
is  used.  This  can  lead  to  discontinuities  in  the  determinant  function  which  precludes  the  use  of  gradient  based  optimization 
algorithms.  Furthermore,  it  induces  multiple  local  minima  that  have  to  be  sorted  through  after  the  optimization  has  been 
performed.  However,  in  any  crystal  system  in  which  the  stiffness  matrix  is  defined  to  be  positive  definite  (a  condition  that 
should  be  true  in  almost  all  cases),  tbe  eigenvalues  are  guaranteed  to  be  complex  and  in  conjugate  pairs  at  v  =  0  (see 
Lemma  1.1  in  [5]).  Furthermore,  there  is  a  velocity,  v;  >  0,  after  which  two  or  more  of  the  eigenvalues  become  real.  Finding 
this  velocity  before  starting  the  iterative  procedure  greatly  enhances  the  stability  of  modern  computational  routines  for 
finding  the  roots  of  a  function.  Thus,  to  begin  the  algorithm  in  an  arbitrary  crystallographic  direction,  the  limiting  velocity 
is  found  first. 

In  [5],  a  method  of  finding  the  limiting  velocity  analytically  for  certain  crystallographic  orientations  was  given.  This 
method  is  also  convenient  to  implement  numerically  so  that  the  limiting  velocity  can  be  found  in  any  arbitrary  direction 
in  any  crystallographic  plane  of  an  arbitrary  crystal.  The  method  starts  by  defining  an  angle,  tp,  used  to  define  rotated  unit 
vectors: 

mi((p)  =  m,  cos{4>)  -F  n,-  sin(0) 
niicp)  =  mi  sin((/))  -F  n,-  cosicp) 

The  vectors  fh  and  n  are  the  m  and  n  vectors,  rotated  by  an  angle  tp  about  the  vector  perpendicular  to  both  of  them,  or 
m  X  n.  This  basically  rotates  the  plane  and  direction  of  propagation  about  the  vector  pointing  out  of  the  paper.  A  rotated 
0,  matrix  is  then  given  as: 


Qik(0)  =  (Cijkt  -  pv^mjmiSikJrhjmi 


(13) 


This  can  be  expanded  by  substituting  (12)  and  written  as: 

Qi/((0)  =  Cijki  (mjmtcos^icp)  -F  (mjn;  -F  m/nj)  cos(</))  sin(0)  -F  njU;  sin^(</))^  -  pv^SikCos^((p) 


(14) 
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Simplifications  were  made  in  this  expression  by  recognizing  that  m,n,-  =  m  ■  n  =  0  and  that  minii  =  n,ni  =  1.  This  expression 
can  then  be  written  in  equation  form  by  using  (3).  This  is  written  as: 

Q(0)  =  D(,/.)-pv2cos2(0)l3>,3 

D(0)  =  (lcos2(0)+  (^R  +  R^^cos(0)sin(</))+Tsin2(</)) 

The  matrix  D(0)  is  known  as  the  acoustical  tensor.  It  is  shown  in  Lemma  1.1  of  [5]  that  this  matrix  is  positive  definite 
for  all  </),  and  thus  has  3  real,  positive  eigenvalues.  It  can  be  seen  from  (15)  that  the  eigenvalues  of  the  matrix  Q  are  then 
expressed  as: 

&■(</>)  =  ^i  (</>)- PV^cos^C^)  (16) 

where  is  the  ith  eigenvalue  of  D((/>).  Lastly,  the  definition  of  the  limiting  velocity  is  the  smallest  velocity  for  which 
the  matrix  Q((/))  becomes  singular  for  some  angle,  (/>.  Furthermore,  Q(</))  is  positive  definite  at  velocities  below  the  limiting 
velocity  but  greater  than  or  equal  to  zero.  Thus,  the  smallest  velocity  for  which  Q(0)  becomes  singular  is  the  point  at  which 
the  smallest  eigenvalue  of  Q((/))  is  zero.  This  problem  can  be  written  as: 


Vi  =  min 

.^€{-90, 90) 


^i(0)  1 

P  |COS(0)| 


(17) 


This  is  a  bounded  ID  optimization  that  requires  only  one  3x3  eigenvalue  solution  per  evaluation  point.  This  problem  can 
be  solved  by  using  a  standard  golden  section  search  or  any  other  found  in  e.g.  [17].  On  a  laptop  PC  this  computation  takes 
roughly  1.7  ms,  which  for  comparison  is  less  than  an  order  of  magnitude  higher  than  the  computation  time  required  for 
one  addition  in  Matlab. 


3.2.  RSW  velocity  calculation 


Once  the  limiting  velocity  is  known,  one  method  of  finding  the  Rayleigh  wave  velocity  is  to  iterate  over  v  until  a  root 
of  the  determinant  is  found.  This  works  well  in  certain  crystal  systems  where  the  determinant  is  well  behaved,  such  as  in 
cubic  crystal  systems.  For  instance,  consider  a  nickel  alloy  material  where  the  stiffness  matrix  is  given  by: 
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The  density  is  assumed  to  be  p  =  8910  kg/m^.  We  seek  the  RSW  velocity  for  a  wave  traveling  in  the  xi -direction  along  the 
crystallographic  plane  given  by  the  Euler  angles  [15  45  8]  in  Bunge’s  active  notation.  This  corresponds  to: 


m  =  {0.9311  0.3514  0.0984)’' 

,  (19) 

n  =  (0.1830  -0.6830  0.7071}’^ 

The  RSW  velocity  will  be  the  velocity  for  which  |B|  x  |B|*  =0.  Plotting  this  quantity  with  respect  to  the  velocity  gives  the 
curve  shown  in  Fig.  1 .  This  curve  was  plotted  from  0  velocity  to  the  limiting  velocity  in  this  direction.  The  RSW  velocity  is 
the  point  at  which  the  function  value  is  zero,  which  is  also  the  global  minimum  of  the  curve  in  the  interval  given  due  to 
taking  the  norm  squared  of  the  determinant  value.  This  minimum  can  be  found  by  using  golden  section  search. 

This  method  works  well  for  cubic  crystal  systems  because  the  determinant  equations  typically  result  in  curves  that  look 
similar  to  Fig.  1.  However,  there  are  pathological  examples  that  make  this  algorithm  unstable  for  certain  directions  in  some 
crystal  systems.  For  instance,  consider  the  hexagonal  Zinc  system  with  stiffness  matrix: 
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(20) 


and  density  p  =  7180  kg/m^.  This  system  gives  the  algorithm  troubles  when  trying  to  compute  RSW  velocity  in  all  directions 
on  the  (010)  plane,  which  corresponds  to  Euler  angles  [0  90  0],  An  example  curve  can  be  seen  in  Fig.  2.  There  are  several 
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Fig.  1.  A  plot  of  the  determinant  of  B  times  the  conjugate  of  the  determinant  of  B  with  respect  to  the  velocity  in  a  cubic  crystal. 


Velocity  (m/s) 

Fig.  2.  A  plot  of  the  determinant  of  B  times  the  conjugate  of  the  determinant  of  B  with  respect  to  the  velocity  in  a  hexagonal  crystal. 


different  minima  in  this  curve  which  seek  to  confuse  the  algorithm.  The  last  minimum  is  the  vaiue  that  corresponds  to  the 
RSW  veiocity  in  this  materiai.  The  curve  couid  be  sampled  many  times  to  find  ail  the  minimum  values,  but  this  leads  to 
computational  inefficiencies  that  are  exacerbated  when  used  within  an  inversion  algorithm.  Thus,  it  is  desirable  to  find  a 
new  method  for  locating  this  minima. 

In  [18],  a  result  is  given  that  is  capable  of  locating  the  RSW  veiocity  without  reiying  on  the  boundary  vaiue  determi¬ 
nant.  In  this  paper,  it  is  noted  that  the  surface  impedance  matrix,  Z,  is  positive  definite  for  0  <  v  <  vr.  This  impiies  that 
the  eigenvalues  of  Z  are  all  positive  untii  the  veiocity  reaches  vr,  at  which  point  one  of  the  eigenvalues  becomes  zero. 
Furthermore,  as  observed  in  [19],  there  is  oniy  one  zero  eigenvalue  at  the  Rayleigh  wave  veiocity,  and  the  eigenvalues  are 
monotonicaiiy  decreasing  as  a  function  of  v  in  the  range  0  <  v  <  vr.  Thus,  if  the  eigenvaiues  are  sorted  in  ascending  order, 
the  first  eigenvaiue  will  have  a  sign  change  at  v  =  vr.  Again,  this  ieads  to  a  quick  and  accurate  method  for  obtaining  the 
veiocity  when  the  limiting  veiocity  is  aiready  known.  The  drawback  to  this  is  that  it  takes  longer  to  compute  the  eigenvalues 
of  2  than  simpiy  using  the  b  vectors  to  compute  the  determinant  of  B.  In  fact,  using  the  standard  algorithms  in  Matlab’s 
fminbnd  function,  the  RSW  veiocity  takes  ~7  ms  to  compute  by  using  the  boundary  value  determinant,  whereas  using  the 
fzero  aigorithm  on  the  smaiiest  eigenvaiue  of  Z  takes  ~20  ms.  Stiii,  this  leads  to  one  forward  simulation  of  180  degrees  of 
orientation  taking  ~3.6  seconds,  which  is  feasible  for  use  in  an  inversion  algorithm. 

4.  Alternate  modes 

In  many  crystals,  more  than  one  surface  wave  mode  may  exist.  This  computationai  routine  is  intended  to  solve  for  elastic 
constants  given  observed  surface  wave  speeds  in  multipie  directions  on  a  single  crystal.  Thus  it  is  important  to  be  abie  to 
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predict  velocities  of  waves  that  are  not  fundamentally  Raylelgh-like  In  nature.  Two  such  waves  will  he  discussed  In  this 
section  corresponding  to  the  crystals  shown  above. 

4.1.  Pseudo-surface  waves  in  cubic  materials 


It  has  been  long  known  that  in  certain  directions  on  certain  surfaces  of  cubic  crystals,  pseudo-surface  waves  will  propa¬ 
gate  that  will  dominate  the  response  over  RSW’s  [20].  A  pseudo-surface  wave  is  one  which  propagates  on  the  surface  for  a 
short  distance.  Portions  of  the  energy  propagate  into  the  sample,  leading  to  decay  in  the  propagation  direction.  Mathemati¬ 
cally,  this  can  be  expressed  as: 


y  _  ag-iklgm-x+pn-x-vt)  ^21) 

where  q  =  1  -I-  hi  is  a  complex  constant  that  Indicates  decay  in  the  propagation  direction.  The  decay  term,  £>  <  0,  is  now  a 
new  constant  that  needs  to  be  estimated  alongside  the  surface  wave  velocity.  Carrying  this  term  through  the  derivations  of 
Stroh’s  eigenvalue  problem  gives  a  new  N  matrix  given  by 


,_r  qT-’R^  T-1 

[-Q'-hq^RT^iR^  -qRT^i 

where 


(22) 


Q' 


3 

ij=t 


(23) 


The  eigenvalues  and  eigenvectors  of  N'  are  used  in  the  same  way  to  predict  the  zeros  of  the  boundary  value  determinant. 
However,  with  the  addition  of  the  complex  terms  in  the  matrix,  the  eigenvalues  no  longer  appear  In  complex  conjugate 
pairs.  Furthermore,  the  third  eigenvalue  will  actually  have  lm(p}  >  0,  which  would  indicate  that  energy  increases  as  depth 
approaches  infinity.  This  increase  is  offset  by  the  decay  in  the  direction  of  propagation. 

When  the  algorithm  encounters  a  condition  where  the  RSW  velocity  is  near  the  limiting  velocity: 


Vr  >  0.99vi  (24) 

it  searches  for  supersonic  solutions  to  the  boundary  value  problem  that  are  just  above  the  limiting  velocity.  The  boundaries 
of  the  optimization  are  moved  to(V(,<v<=vjv).  and  optimization  performed  again.  In  the  current  iteration  of  the  algo¬ 
rithm,  vn  is  chosen  to  be  the  point  at  which  there  are  more  than  2  real  eigenvalues,  which  can  be  found  by  an  Iterative 
search  routine. 

The  solutions  to  this  new  optimization  problem  may  not  be  zeros,  but  they  are  minima  of  the  boundary  value  deter¬ 
minant.  If  they  are  not  zeros,  the  decay  term  is  added  to  the  eigenvalue  problem  and  the  search  is  continued  on  both 
b  and  v.  The  advantage  to  this  approach  is  that  this  second  optimization  is  performed  only  in  the  directions  at  which 
a  pseudo-surface  wave  is  likely  to  exist.  As  the  algorithm  progresses  through  more  directions,  the  previously  calculated 
pseudo-surface  wave  velocity  can  be  used  as  an  initial  guess  to  make  the  second  optimization  converge  very  quickly. 


4.2.  Alternate  surface  wave  modes  in  non-cubic  materials 


Surface  waves  propagating  on  the  basal  plane  of  hexagonal  crystals  do  not  behave  like  typical  Rayleigh  waves.  The 
dominant  surface  wave  mode  on  this  specific  crystallographic  plane  Is  a  supersonic,  2-component  wave  (according  to  the 
classification  system  from  [21]).  The  velocity  of  this  wave  was  predicted  analytically  by  [22[.  As  the  cut-plane  tilts  away  from 
the  basal  plane,  this  supersonic  mode  eventually  gives  way  to  the  traditional  subsonic  RSW  with  all  three  components.  In 
the  configuration  in  which  the  free  surface  normal  lies  in  the  Basal  plane,  the  surface  waves  that  propagate  in  any  direction 
are  purely  traditional  Rayleigh  waves  and  their  velocities  have  been  predicted  in  certain  directions  by  many  authors  (e.g. 
[23]).  The  condition  for  the  existence  of  this  alternate  propagation  mode  in  hexagonal  crystals  Is  the  same  as  that  for 
the  PSW  In  cubic  crystals,  and  the  method  of  finding  the  velocity  Is  also  the  same.  The  calculations  were  verified  on  a 
hexagonal  crystal  analyzed  previously  In  [24[.  In  that  work,  the  authors  performed  computations  for  a  number  of  different 
crystals  with  the  reference  plane  of  the  surface  wave  coincident  with  a  plane  of  symmetry  of  a  the  crystal  system.  Their 
plots  show  propagation  of  surface  waves  along  a  surface  with  the  angle  between  the  surface  normal  and  the  main  axis,  or 
c-axls,  of  the  hexagonal  material  varying  between  0  and  90  degrees.  The  algorithms  from  the  current  work  were  used  to 
calculate  this  curve  for  cadmium  telluride,  and  the  data  is  shown  in  Fig.  3.  This  plot  shows  the  same  surface  wave  velocity 
as  “Fig.  11”  from  [24]. 

To  test  this  algorithm  on  more  complex  material  systems,  a  monoclinic  material  was  analyzed  as  well.  In  [25],  the 
authors  discussed  propagation  of  surface  waves  when  the  plane  of  symmetry  in  monoclinic  materials  was  spanned  by  the 
normal  to  the  free  surface  as  well  as  the  direction  of  propagation,  n  and  m,  in  the  notation  used  in  this  paper.  They  gave 
an  in-depth  analysis  of  the  theory  of  wave  propagation  in  this  special  case  and  showed  a  number  of  numerical  examples. 
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Angle  Between  Surface  Normal  and  c-Axis  (degrees) 

Fig.  3.  Plot  of  the  surface  wave  velocity,  the  limiting  velocity,  and  the  fast  and  slow  shear  wave  velocities  for  cadmium  telluride  with  the  c-axis  of  the 
hexagonal  crystal  in  the  plane  spanned  by  m  and  n.  The  x-axis  is  the  angle  of  the  primary  axis  of  the  hexagonal  crystal  with  respect  to  the  normal  of  the 
surface. 


Fig.  4.  Plot  of  the  surface  wave  velocity,  the  limiting  velocity,  and  the  fast  and  slow  shear  wave  velocities  for  aegirite-augite  with  the  surface  waves 
propagating  in  the  plane  of  symmetry. 


For  the  purpose  of  this  work,  one  crystal  that  they  analyzed,  aegirite-augite,  was  used  to  calculate  the  surface  wave  velocity 
and  compare  with  the  results  shown  in  this  previous  work.  The  curves  shown  in  that  work  were  calculated  by  keeping  n 
and  m  in  the  plane  of  symmetry  and  rotating  them  by  an  angle,  a,  about  the  vector  normal  to  them,  essentially  changing 
the  angle  of  inclination  of  the  stress  free  boundary  relative  to  the  crystal  system.  The  results  of  the  computation  are  shown 
in  Fig.  4.  Once  again,  the  predicted  surface  wave  velocity  is  the  exact  same  curve  as  that  shown  in  the  paper  by  Chadwick 
and  Wilson. 

5.  Verification  of  the  numerical  scheme 

The  algorithms  for  calculating  surface  wave  velocities  discussed  in  the  previous  section  were  verified  using  other  numer¬ 
ical  methods.  This  section  shows  these  results  as  well  as  verification  in  the  context  of  inverse  algorithms  for  calculation  of 
the  elastic  constants. 

5.1.  Calculation  of  surface  wave  velocity  on  certain  cut-planes 

A  large  set  of  experimental  data  sets  on  multiple  different  cut  planes  in  a  cubic  nickel  crystal  was  given  in  [2],  In  this 
paper,  two  specific  planes  from  [2]  were  picked  to  determine  the  agreement  with  experimental  data.  In  their  paper,  Li  et 
al.  gave  the  crystallographic  planes  as  Miller  indices,  so  the  indices  were  converted  to  Euler  angles  for  use  in  the  code. 
The  same  material  properties  as  those  found  in  the  paper  were  used  to  generate  the  curves  (cn  =  235  GPa,  cu  =  L42  GPa, 
cn  =  131  GPa,  and  p  =  8720  kg/m^).  The  velocities  generated  by  the  code  discussed  previously  are  given  in  Fig.  5.  The 
first  plot.  Fig.  5a  shows  the  velocities  on  the  plane  (0.28  0.1  1)  and  the  second  is  for  the  plane  (0  0.25  1).  While  the  tilt  in 
plane  is  different,  the  calculated  velocities  are  the  same  as  those  shown  in  the  paper,  which  also  matched  the  experiments 
that  were  performed  in  that  work.  These  figures  take  an  average  of  ~3  s  to  generate  on  a  laptop  computer  with  an  Intel 
Core™i7-3630QM  processor  and  16  GB  of  DDR3  RAM. 
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Rayleigh  Wave  Velocity  on  Surface  Rayleigli  Wave  Velocity  on  Surface 


(a)  (b) 

Fig.  5.  Polar  plots  showing  the  RSW  velocity  and  the  PSW  velocity  on  the  nickel  crystal  used  in  [2].  (a)  Shows  the  velocity  on  the  (0  0.25  1)  plane  and 
(b)  shows  the  velocity  on  the  (0.28  0.1  1)  plane. 


Rayleigh  Wave  Velocity  on  Surface 


Fig.  6.  Plot  of  the  wave  velocities  on  the  plane  used  in  the  inverse  problem  verification. 


5.2.  Verification  of  the  inverse  problem  on  cubic  materials 


The  inverse  problem  was  tested  with  simulated  data  to  verify  convergence  and  test  the  sensitivity  to  the  initial  guess 
of  the  elastic  constants.  Velocity  curves  were  generated  using  the  elastic  constants  given  in  §3.2  on  the  cubic  crystal  given 
there  for  the  plane  described  by  the  Euler  angles  [15  45  8].  Noise  was  added  to  the  data  and  the  inverse  problem  was 
run  for  multiple  initial  guesses  of  elastic  constants  to  determine  the  convergence  of  the  algorithm  with  increasingly  poor 
starting  points.  The  velocity  profile  at  this  orientation  is  shown  in  Fig.  6. 

In  an  experiment,  the  velocity  of  sound  is  typically  not  measured  directly.  Instead,  often  the  arrival  time  of  a  wave  is 
measured  at  one  point  on  a  sample  surface  and  then  again  on  a  different  point  separated  by  a  distance.  Ad.  Thus,  the 
noise  in  the  experiment  will  likely  be  on  the  measured  arrival  time.  In  lieu  of  this,  synthetic  data  was  generated  by  adding 
independent  and  identically  distributed  (IID)  noise  to  the  calculated  change  in  arrival  time  for  the  waves  in  each  direction: 


synth  _  calc  ,  p.  -  Af  _l  c. 

t.  -  tj  +  e,  -  +  e, 


e,  ~  Normal(0,  a^) 


(25) 


The  calculated  arrival  time,  is  determined  by  calculating  the  velocity  in  the  ith  direction,  and  dividing  an  assumed 
spacing  between  measurement  points  by  this  velocity.  The  assumption  that  the  noise  in  the  data  is  IID  and  normal  is 
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Simulated  Data  with  High  Noise  Level 


Fig.  7.  The  arrival  times  as  a  function  of  angle  away  from  the  sample  X  axis,  with  the  solid  line  showing  the  computed  arrival  times  and  the  scatter  plot 
showing  the  arrival  times  with  noise  added. 


assumed  to  be  valid  in  the  experiment,  though  possible  deviations  should  be  addressed  if  the  assumptions  do  not  appear 
valid  from  analysis  of  the  data.  This  allows  the  algorithm  to  take  advantage  of  the  typical  nonlinear  least  squares  estimates: 


argmin 


m=l 


uSynth 


subject  to  c ' 


lower 

O' 


<  Cij  i 


Ad  Y 

Vm(Cij)J 

upper 


(26) 


The  change  in  arrival  time  for  the  plane  used  in  this  section  at  a  separation  distance  of  .01  m  is  given  in  Fig.  7.  The 
simulated  data  is  plotted  on  top  of  the  calculated  arrival  times.  In  this  study,  a  standard  deviation  of  3  x  10“®  s  was  used 
to  generate  the  noise  in  the  data.  This  number  was  chosen  based  on  previous  experience  with  peak  fitting  for  finding  the 
time  of  arrival  of  Rayleigh  surface  waves.  The  experimental  setup  in  the  lab  was  discussed  in  [26].  Though  no  number  was 
given  there  for  the  standard  deviation  of  time  of  arrival,  3  x  10“®  s  has  been  typical  for  this  experiment. 

Two  aspects  of  the  forward  problem  make  the  inverse  problem  more  challenging.  Scaling  is  important  factor  in  any 
optimization  problem.  For  instance,  the  results  of  the  simulation  are  in  seconds,  but  the  time  of  flight  of  waves  in  an 
ultrasound  experiment  are  generally  on  the  order  of  ps.  Therefore,  the  objective  function  shown  in  (26)  is  multiplied  by  a 
factor  of  le6.  Furthermore,  elastic  constants  are  generally  on  the  order  of  lell,  so  the  parameters  of  the  inverse  function  are 
divided  by  this  factor.  This  ensures  that  the  scales  of  the  objective  function  and  parameters  are  similar.  Another  confounding 
factor  is  that  the  gradients  of  the  problem  are  not  currently  calculated  analytically.  Furthermore,  the  objective  function 
is  itself  a  minimization  problem  which  is  subject  to  non-negligible  noise  due  to  the  convergence  criteria  of  the  inner 
optimization  problem.  Finite  difference  methods  for  gradient  computation  tend  to  enhance  the  strength  of  this  noise  to 
the  point  that  derivative  calculation  with  finite  difference  is  very  unreliable.  A  combination  of  finite  difference  for  the 
eigenvalue  problem  and  analytical  derivatives  for  the  optimization/root  finding  algorithm  can  be  used  to  reliably  compute 
the  sensitivities,  but  this  is  at  considerable  computational  cost.  Thus,  for  this  work,  a  non-gradient  based  Nelder-Mead 
optimization  algorithm  was  used  to  perform  the  optimization  in  the  inverse  problem. 

Another  important  observation  is  that  the  inverse  problem  is  generally  more  well  behaved  for  cubic  crystals  when  the 
parameter  estimates  are  performed  on  the  alternative  but  equivalent  parameterization: 


(cii  -F  2ci2)(cii 
Cll  +  Ci2 


V  = 

a  = 


Cl2 

Cll  +  Ci2 
2C44 

Cll  —  Ci2 


Cu) 


(27) 


This  parameterization  is  essentially  the  Young’s  modulus,  Poisson  ratio,  and  the  Zener  anisotropy  ratio  of  the  material, 
respectively.  The  sensitivity  of  the  velocity  with  respect  to  these  two  parameterizations  were  computed  at  the  known  elastic 
constants  and  are  plotted  in  Fig.  8.  It  can  be  seen  in  these  plots  that  the  sensitivity  to  cn  and  C]2  are  both  very  similar 
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Sensitivity  of  Dominant  Modes  Sensitivity  of  Dominant  Modes 


(a)  (b) 


Fig.  8.  Figures  showing  the  sensitivity  of  the  dominant  wave  mode  velocity  with  respect  to  (a)  the  cy  parameterization  and  (b)  the  E.  v,  and  a  parameter¬ 
ization. 


Sensitivity  of  Dominant  Modes 


Fig.  9.  The  sensitivity  of  the  dominant  wave  mode  velocity  with  respect  to  the  eigenvalue  parameterization. 


for  all  the  directions  of  propagation,  which  would  tend  to  indicate  identifiability  issues  with  these  two  parameters.  This  is 
not  an  issue  with  the  £,  v,  and  a  parameterization.  An  alternative  parameterization  could  be  to  use  the  eigenvalues  of  the 
stiffness  matrix  as  the  parameters  for  the  inverse  problem,  or: 

Al  =  Cii  —  Ci2 

A2  =  Cii-|-2ci2  (28) 

Ag  =  2C44 

This  parameterization  has  the  added  advantage  of  simply  being  able  to  place  bounded  constraints  on  the  optimization  to 
keep  the  eigenvalues  positive,  thus  enforcing  positive  definite  stiffness  matrices.  However,  the  sensitivities  (shown  in  Fig.  9) 
have  similar  issues  as  those  from  the  usual  Cjj  parameterization  between  Ag  and  >.3.  Thus,  the  £,  y,  and  a  parameterization 
is  used  in  the  inverse  problem,  and  the  results  of  the  inverse  are  converted  hack  to  the  usual  Cjj  after  an  optimal  solution 
has  been  found. 

The  inverse  problem  was  run  with  multiple  different  initial  design  points.  The  known  values  of  the  elastic  constants  that 
generated  the  data  were  [2.41, 1.47, 1.26]  x  10^’  Pa.  Initial  design  points  were  selected  randomly  by  assuming  a  uniform 
distribution  in  a  reasonable  range  of  £,  v  and  a.  The  table  of  the  initial  design  points  is  shown  in  Table  1.  The  final  values 
for  each  starting  point  are  shown  in  Table  2.  It  is  clear  that  the  initial  designs  all  underestimated  the  values  of  cy,  the 
most  egregious  of  which  are  points  4-7.  However,  the  only  design  points  that  did  not  converge  to  the  correct  minimum  are 
points  2,  3,  and  8.  These  points  had  initial  values  that  were  relatively  close  to  the  true  values,  with  3  being  the  closest  out 
of  all  the  initial  designs.  This  strange  behavior  is  likely  due  to  the  algorithm  used  for  the  optimization  not  being  able  to 
effectively  linearize  the  problem  around  these  design  points.  This  is  because  of  the  numerical  error  in  the  inner  optimization 
loop  from  Equation  (10).  However,  the  values  of  the  objective  function  are  three  orders  of  magnitude  lower  for  points  1, 
4-7,  and  9-10  than  they  are  for  the  troubled  values,  and  it  is  clear  that  the  algorithm  can  actually  be  used  effectively  in  the 
inverse  problem  provided  multiple  initial  designs  are  sampled.  Strides  should  be  taken  to  improve  the  gradient  estimation 
so  that  a  gradient-based  approach  for  the  inverse  problem  can  be  used  rather  than  the  linearized  simplex  algorithms. 
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Table  1 

Starting  points  attempted  for  inverse  (lell  Pa). 


Point  number 

cu 

Cu 

C44 

1 

1.9201 

1.0669 

0.9903 

2 

2.0366 

0.9923 

1.1760 

3 

2.3138 

1.3255 

1.1245 

4 

1.5680 

0.7683 

1.0896 

5 

1.5746 

0.7735 

1.0384 

6 

1.5329 

0.7445 

0.9812 

7 

1.4871 

0.7323 

0.8413 

8 

1.8068 

0.8971 

1.0227 

9 

1.9947 

0.9992 

1.0723 

10 

2.1144 

1.1001 

1.3895 

Table  2 

Final  values  for  each  starting  point  (lell  Pa). 

Point  number 

Cll 

Cu 

C44 

1 

2.3365 

1.3942 

1.2683 

2 

4.2637 

3.3114 

1.1234 

3 

4.1307 

3.1769 

1.1174 

4 

2.3365 

1.3941 

1.2683 

5 

2.3364 

1.3940 

1.2683 

6 

2.3399 

1.3976 

1.2681 

7 

2.3365 

1.3941 

1.2683 

8 

4.8179 

3.8537 

1.0951 

9 

2.3364 

1.3940 

1.2683 

10 

2.3365 

1.3941 

1.2683 

6.  Summary 

A  new  numerical  scheme  was  developed  for  calculating  the  velocity  of  the  dominant  surface  wave  mode  in  anisotropic 
elastic  media.  The  algorithm  was  applied  to  hexagonal  and  cubic  crystals  and  found  to  be  accurate.  The  method  takes  3 
seconds  to  build  an  entire  velocity  cui've  for  181  angles  of  wave  propagation  on  an  arbitrary  face  of  a  crystal,  which  is  a 
significant  improvement  in  the  amount  of  time  previous  implementations  of  numerical  schemes  have  taken.  The  algorithm 
was  applied  in  the  numerical  inverse  problem  for  calculating  the  elastic  constants  of  a  cubic  crystal  with  data  generated 
by  adding  IID  noise  to  the  time-of-arrivals.  The  algorithm  converged  to  the  correct  minimum  in  the  inverse  problem  for 
seven  out  of  10  initial  designs,  and  improvements  in  the  gradient  computation  could  result  in  the  use  of  gradient-based 
optimization  algorithms  which  do  not  have  as  much  sensitivity  to  the  initial  design  point.  This  remains  as  future  work. 
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